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ABSTRACT 

Context. The Kepler mission’s discovery of a number of circumbinary planets orbiting close (Up <1.1 au) to the stellar binary raises 
questions as to how these planets could have formed given the intense gravitational perturbations the dual stars impart on the disk. The 
gas component of circumbinary protoplanetary disks is perturbed in a similar manner to the solid, planetesimal dominated counterpart, 
although the mechanism by which disk eccentricity originates differs. 

Aims. This is the first work of a series that aims to investigate the conditions for planet formation in circumbinary protoplanetary 
disks. 

Methods. We present a number of hydrodynamical simulations that explore the response of gas disks around two observed binary 
systems: Kepler-16 and Kepler-34. We probe the importance of disk viscosity, aspect-ratio, inner boundary condition, initial surface 
density gradient, and self-gravity on the dynamical evolution of the disk, as well as its quasi-steady-state profile. 

Results. We find there is a strong influence of binary type on the mean disk eccentricity, ej, leading to = 0.02 - 0.08 for Kepler-16 
and grf = 0.10 - 0.15 in Kepler-34. The value of a-viscosity has little influence on the disk, but we find a strong increase in mean disk 
eccentricity with increasing aspect-ratio due to wave propagation effects. The choice of inner boundary condition only has a small 
effect on the surface density and eccentricity of the disk. Our primary finding is that including disk self-gravity has little impact on 
the evolution or final state of the disk for disks with masses less than 12.5 times that of the minimum-mass solar nebula. This finding 
contrasts with the results of self-gravity relevance in circumprimary disks, where its inclusion is found to be an important factor in 
describing the disk evolution. 

Key words, methods: numerical - hydrodynamics - planets and satellites: formation - protoplanetary disks - binaries: close 


1. Introduction 


Extrasolar planets in circumbinary (p-type) orbits around short- 
period stellar binary systems are a recent addition to the increas¬ 
ingly diverse range of planetary characteristics found outside our 
solar system. We now know of nine short-period planet^ with 
semi-major axis Up < 1.1 au, placing them close to the binary 
barycentre and under the influence of perturbative forces from 
the secondary star. 

The stability of circumbinary planet orbits was studied well 
before the first confirmed exoplanet detection in 1992 ( |Wol-| 
szczan & Erail||1992|l by |Heppenheimer ( 1978| l. More recently, 
Holman & Wiegert] ^1999| ) performed a parameter space study 
of the long-term stability of planets in P-type orbits. They calcu¬ 
lated a range of orbital radii that allows for dynamical stability 
for a range of binary eccentricities and semi-major axes. Inter- 


‘ Kepler-16(AB)b ( |Doyle et al.||201l| l , Kepler-34(AB)b & 35(AB)b 
( [Welsh et al.||2012|l, Kep ler-38(AB)b ( jOrosz et al.||2012bj, Kepler- 
47(AB)b,c (jOrosz et al.||2012aj, Kepler-64(AB)b JSchwamb et al. 
IoT 3)|, Kep ler-413(AB)b ( [Kostov et al.|2014| l and KIC 9632895 ( [Welsh 
et al. 2014b. 


estingly, with the exception of Kepler-47(AB)c, all known cir¬ 
cumbinary planets lie just outside their inner most stable orbit, 
Uc. One plausible explanation for this occurrence is that at some 
point in their evolution, these planets underwent a period of mi¬ 
gration that was halted. 


[Pierens & Nelson] ( 2008| l showed that Saturn- and Jupiter- 
mass planets undergo different interactions with the binary that 
can result in either stabilisation of the planetary orbits or its re¬ 
moval from the system via ejection or scattering. Specifically 
they find that after a period of inwards migration, Saturn-mass 
planets have a stable evolution. A torque reversal, caused by an 
eccentricity increase from interaction with the binary, leads to 
stable outwards migration. However, Jupiter-mass planets often 
become trapped in 4; 1 mean motion resonance with the binary, 
which causes significant increases in planetary eccentricity ( |Nel-[ 
son|2003| l. Eventually, close encounters with the stars at periap- 
sis lead to outward scattering or complete ejection. The latter 
effect may explain why all Kepler circumbinary planets discov¬ 
ered so far are sub-Jupiter mass, a result difficult to explain by 
observational biases since larger planets are easier to detect. 
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Fig. 1. Surface density maps, 1,20, for non self-gravitating models around Kepler-16 (top row) and Kepler-34 (bottom row) at 16,000 Pab- Frames 
a-d and e-h correspond to simulations A-D and E-H in Table[T] respectively. 


Although all the observed circumbinary planets are currently 
in stable orbits, it is not obvious that they formed in situ. Several 
studies have been carried out to try to understand how the pro¬ 
cesses behind planet formation and evolution can occur in spite 
of an intense gravitational influence on the protoplanetary disk 
from the central binary. One of the key hurdles in forming solid 
cores and terrestrial planets in circumbinary disks is explain¬ 
ing how rocky, metre- to kilometre-sized planetesimals undergo 
collisions that result in growth rather than erosion. It has been 
shown in previous work that eccentricity forcing from the binary 
on the protoplanetary disk can drive up relative velocities which 
not only inhibits accretion but can lead to energetic and highly 
erosive impacts. In particular, [Lines et’aL] ( |2014| ) showed that 
even massive planetesimals are often disrupted despite having a 
large gravitational binding energy. Their work corroborates that 
of jPaardekooper et al.| ( |2012l l and |Meschiari| ( |2012a|b| l suggest¬ 
ing that sustained planetesimal accretion is unlikely and thus, 
the observed Kepler circumbinary planets did not form in situ. 
The most likely explanation for the presence of these planets is 
the migration of planetary cores that form at larger orbital radii 
where the gravitational forcing from the binary is of significantly 
lower magnitude. 


A key component missing from many of the A-body studies 
of planetesimal evolution in circumbinary disks is the inclusion 
of a gas disk which is typically 100 times the mass of its solid 
counterpart. As is the case with the rocky planetesimals, the gas 
is also perturbed by the binary, raising the eccentricity of the 
disk. The eccentricity in the gas disk is established through para¬ 
metric instabilities originating from non-linear mode coupling 


between the eccentric m - 
centricity and the binary pote 

1 modes of both the initial disk ec- 

mtial (Papaloizou 

2002[[Hayasaki & 

Okazaki|2009 [Lubow|1991 

and well as the c 

irect driving from 

the binary m = 1 potential 

Lubow & Artymowicz||2000|l. This 
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coupling produces an m = 2 spiral wave at the 3:1 Lindblad res¬ 
onance which increases disk eccentricity via the outwards trans¬ 
port of angular momentum. The overall result is an eccentric, 
asymmetric disk precessing about the binary (jKley & Haghigh 


ipour| 


arzari et ^2013| Pelupessy & Portegies Zwart 


2013[ |Pierens & Nelson|2Cil3| (hereafter referred to as PN13)). 


The asymmetric gas disk interacts with the solid planetesimal 
disk through gas drag and the time dependent gravitational po¬ 
tential of the gas disk. The latter effect is particularly important 
when large asymmetries in the gas surface density arise (see Fig¬ 
ure [^. 

The role of gas drag has been probed in several studies which 
found that differential orbital phasing caused relative velocities 
to remain small between planetesimals of comparative mass, 
thus, orbit crossing will occur between planetesimals of differ¬ 
ent sizes ( [Scholl et al.|2007j|Marzari et al.|2008l l. The influence 
of disk self-gravity on the evolution of the system is less well 


explored. Recent studies such as PN13 and Kley & Haghigh- 


[ipour[ ( |2014[ i choose not to include self-gravity in their Simula 
tions. The former justifies this approximation by confirming the 
Toomre parameter. 


Q 


Cs^ 

ttGI’ 


( 1 ) 


where Cs is the sound speed, Q is the angular velocity, G is 
the gravitational constant, and 2 is the surface density, satis¬ 
fies the Toomre stability criterion {Q > 1) at all disk radii for 
their chosen disk mass. However, even in self-gravitating disks 
of relatively low mass, low-frequency global modes exist that do 
not have a counterpart in non self-gravitating disks (Papaloizou 
2002| l. Moreover, the study by [Marzari et al. ( 2009| l shows that 
self-gravity in circumprimary disks can be important even in 
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Simulation 

Binary System 

h 

a 

Fo {Malau^) 

Self-Gravity 

Boundary 

CTl 

A 

Kepler-16 

0.03 

10-" 

1.0 X 10-"^ 

No 

Rigid 

1.5 

B 

Kepler-16 

0.05 

10-3 

1.0 X 10-"^ 

No 

Rigid 

1.5 

C 

Kepler-16 

0.03 

10-2 

1.0 X 10-"^ 

No 

Rigid 

1.5 

D 

Kepler-16 

0.05 

10-2 

1.0 X 10-"^ 

No 

Rigid 

1.5 

E 

Kepler-34 

0.03 

10-3 

1.0 X 10-4 

No 

Rigid 

1.5 

F 

Kepler-34 

0.05 

10-3 

1.0 X 10-4 

No 

Rigid 

1.5 

G 

Kepler-34 

0.03 

10-2 

1.0 X 10-4 

No 

Rigid 

1.5 

H 

Kepler-34 

0.05 

10-2 

1.0 X 10-4 

No 

Rigid 

1.5 

I 

Kepler-16 

0.03 

10-2 

1.0 X 10-4 

Yes 

Rigid 

1.5 

J 

Kepler-16 

0.05 

10-3 

1.0 X 10-4 

Yes 

Rigid 

1.5 

K 

Kepler-34 

0.03 

10-3 

1.0 X 10-4 

Yes 

Rigid 

1.5 

L 

Kepler-34 

0.05 

10-3 

1.0 X 10-4 

Yes 

Rigid 

1.5 

M 

Kepler-16 

0.05 

10-2 

1.0 X 10-4 

No 

Open 

1.5 

N 

Kepler-34 

0.05 

10-3 

1.0 X 10-4 

No 

Open 

1.5 

O 

Kepler-16 

0.05 

10-2 

5.0 X 10-4 

Yes 

Rigid 

1.5 

P 

Kepler-16 

0.05 

10-3 

1.0 X 10-3 

Yes 

Rigid 

1.5 

Q 

Kepler-16 

0.05 

10-3 

2.5 X 10-3 

Yes 

Rigid 

1.5 

R 

Kepler-16 

0.05 

10-3 

5.0 X 10-4 

No 

Rigid 

1.5 

S 

Kepler-16 

0.05 

10-3 

1.0 X 10-3 

No 

Rigid 

1.5 

T 

Kepler-16 

0.05 

10-3 

2.5 X 10-3 

No 

Rigid 

1.5 

U 

Kepler-16 

0.05 

10-2 

1.0 X 10-4 

Yes 

Rigid 

0.5 


Table 1. Parameter setup of each simulation A through U. A-H test the effect of varying aspect ratio, alpha-viscosity and central binary type on 
non self-gravitating disks. I-L test the inclusion of self-gravity at standard surface density. M-N investigate the inclusion of an open boundary. 0-T 
look at the relevance of self-gravity in more detail by increasing the surface density of the disk. U tests the dependency of the disk evolution on 
the surface density profile exponent, a^. 


low-mass disks. Therefore, an investigation into whether or not 
self-gravity is required to fully describe a circumbinary disk, 
even at lower masses, is required. 

In this paper, we aim to explore a number of parameters 
characterising a gaseous circumbinary disk in an attempt to find 
a suitable surface density profile for both the Kepler-16 and 
Kepler-34 systems. The resulting quasi-steady-state surface den¬ 
sity profiles will be used to parameterize the gas disk potential 
for use in our A-body simulations of planetesimal growth. This 
next work will be presented in a second paper, part 11. Section 
[^presents the numerical method and model parameters. In Sec¬ 
tion and we discuss the results and implications of the 21 
simulations of the Kepler-16 and Kepler-34 systems. Finally, we 
summarise our findings in the conclusion in Section]^ 


2. Numerical methods 


The central stellar binary potential is calculated in a fixed steady 
state orbit centred on the binary barycentre, meaning that the 
stars do not respond to the gas disk ( jPierens & Nelson||2007| l. 
The stars, with mass ratio q - MbIMa are first initialised at 
their respective apoapses, such that the Cartesian coordinates are 
initially zero in the Y direction for both stars, Ya & Yb — 0, and 
the X position is 




and 


I + q 


Xb^ 


\ + q 


( 2 ) 


(3) 


where Ma and Mb are the mass of star A and B, respectively, 
and the distance of a star from its focus, r, is determined from 
the shape equation 


r - at, 


1 


1 - et cos 0 I ’ 


(4) 


where at, is the binary semi-major axis, et is the binary eccentric¬ 
ity, and the angle of the star from apoapsis is 0. Using Kepler’s 
2nd law to determine the fractional area swept out at any given 
time, an iterative procedure is invoked to calculate 0(t) and thus, 
the position vectors of each star. The stellar binary parameters 
are given in Table 

To simulate the hydrodynamical evolution of the circumbi¬ 
nary disks we use a modified version of the Eulerian fluid 
code FARGO (|Masset||2000) l, called FARGO-ADSG ( |Baruteau| 
& Masset|2008 1, which includes optional disk self-gravity. The 


hydrodynamical equations are solved on a polar mesh which is 
composed of = 512 equally divided azimuthal sector cells, 
and Hr = 395 logarithmically spaced radial cells. Logarithmic 
spacing in the radial direction is a necessary condition for the 
fast fourier transform algorithm in the self-gravity calculation, 
but is also a preference for achieving the highest resolution clos¬ 
est to the binary ( |Baruteau & Masset|2008[ ). 

The majority of our runs use a rigid, reflecting inner bound¬ 
ary condition which is fixed at 0.345 au with the rigid outer 
boundary at 4.0 au. In a couple of simulations we have used 
standard outflow boundaries at the grid’s inner and outer edges, 
where the zero-gradient condition has been extended to the gas 
azimuthal velocity. This is a necessary implementation since 
there is no well defined equilibrium between the central grav- 
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Kepler-16 


Kepler-34 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

Vr ( au / yr) 

Fig. 2. Radial velocity maps at 16,000 Pab for a) Kepler-16 a = 10“^, h = 0.03, b) Kepler-16 a = 10“^, h = 0.05, c) Kepler-34 a = 10“^, h = 0.03 
and d) Kepler-34 a = 10-^ h = 0.05. 


System 

q 

eb 

Qb (au) 

Kepler-16 

0.29 

0.16 

0.22 

Kepler-34 

1.0 

0.52 

0.23 


Table 2. Stellar binary parameters. 


ity from the binary and the opposing centrifugal force and pres¬ 
sure gradient. In all simulations the inner disk edge is truncated 
by the binary causing a steep positive density gradient which is 
numerically difficult to model smoothly, thus, we employ a gap 
function, from|Gunther et al.|(|2004|l to better model the ini¬ 
tial surface density prohle in the tiiner disk region: 


fs 


gap 


1 + exp 


R-Rs 

0.\R.c 


(5) 


where Rgap - 2.5ah is the estimated size of the gap (Artymowicz 
|& Lubow(l994[ ) and R is the orbital radius. 

The simulated gas disk is simplified by using the isother¬ 
mal disk approximation, thereby minimising the number of un¬ 
knowns. We leave the inclusion of an energy equation to account 
for the thermal evolution of the disk for a future paper. 

The initial surface density of the gas disk follows S(/?) = 
fgap^{)R~°'^ from PN13, where So is the surface density at 1 au, 
R is the radial position in au, and defines the initial surface 
density gradient. For the majority of the simulations (Table 
simulations A-N & U) we assui ne a half minim um-mass solar 
nebula, thus, Sq = 10“"^ Mo/au^ ( Hayashijl981 1 . In simulations 
0-T So is increased by a factor of 5, 10 and 25. The density gra¬ 
dient, aj is set to 1.5 in all cases except the final simulation, U, 
where a shallower value of 0.5 is used for comparison to previ¬ 
ous work ( [Marzari et al.|20l'3] l. In all simulations, the total cen¬ 
tral stellar mass is normalised to 1.0 Mq but the surface density 
of the disk is not scaled to account for the physical mass differ¬ 
ence between Kepler-34 and Kepler-16. Therefore, although the 
total physical mass of the disk is fixed, the disk mass relative to 
combined stellar mass changes. 

To avoid numerical instabilities caused by extremely low 
density fluid cells in the inner cavity due to the torque of the 
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binary on the disk inner edge, a density floor is added such that 
the minimum value, is set to 10“® of the initial density. This 
means that mass is not strictly conserved in regions of very low 
density. 

The aspect ratio, h - R, where H is the pressure scale height, 
is constant across the disk and varies between each model from 
0.03 to 0.05. The aspect ratio does not determine the physical 
thickness of the disk, since the simulations are two dimensional, 
but defines the sound speed, Cs - Vkh, where Vk is the Keplerian 
speed. Turbulence in the disk, a likely contribution of the MRI 
effect (Balbus & Hawley 1991|l , is explored by varying an a- 
viscosity between 10“^ and 10“^. 

Since our disks are low in mass, we expect them all to be 
linearly stable according to the Toomre stability criterion, Q> \ 
( |Toomre|1964| l. For our half-MMSN density models, even near 
the inner edge of the disk where the density is at a maximum, the 
minimum Toomre value, Qmini is 140. However, on consideration 
of the results of |Marzari et al. (20091 that show self-gravity plays 
a significant role in circumprimary disks, self-gravity is included 
in simulations I-L and 0-Q. 

Over 21 simulations, we explore the effect that changing 
the stellar binary parameters, aspect ratio, a-viscosity, boundary 
condition and inclusion of disk self-gravity has on the evolution 
and quasi steady-state (QSS) profile of the disk. All simulations 
are evolved for 16,000 binary orbits, Pab, equivalent to 1,700 
orbits at 1 au, with the exception of M & N which run for 4,000 
Pab to test the effect of the chosen inner boundary condition on 
the initial response of the disk. 


3. Results 

3.1. Non self-gravitating disks 

3.1.1. Kepler-16 

We first look at the response of a non self-gravitating disk to 
a stellar binary with parameters chosen to match the observed 
properties of the system Kepler-16(AB). The azimuthally and 
time averaged quasi-steady state surface density is shown in Fig¬ 
ure Hh- The density maximum or peak, 'Lpeak, indicates the loca¬ 
tion of the truncated inner edge, although the eccentric morphol¬ 
ogy of the disk requires a full 2-dimensional density map, I, 2 d 
(Figure^, to identify the edge precisely as a function of azimuth 
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Kepler-16 


Kepler-34 



Fig. 3. Time-averaged quasi-steady state surface density profiles (a & b), E, and time evolution of disk mean eccentricity (c & d), ej, and time- 
averaged radial eccentricity profiles e & f), e, for Kepler-16 (Models A - D) and Kepler-34 (Models E - H) under non self-gravitating conditions. 
Time-averaging is over the last 1,000 binary orbits of the simulations. 
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angle. Therefore, we will consider the location of 'Zpeak to be at 
the average truncation radius, r,. The value of r, lies between 1.0 
and 1.4 au which is consistent with what was found by PN13. 

Increasing the aspect ratio from 0.03 to 0.05 raises the disk 
eccentricity (Figures &[^), particularly for a - 10“^, which 
can also be seen in the surface density profiles by an increase 
in the value of f,. The more circular form of the disk in model 
A leads to a higher averaged surface density, since more mass 
orbits at a similar orbital radius. 


The eccentric nature of the disk is better shown by the time 
evolution of the mean disk eccentricity (Figure]^), dd, which we 
define according to Pierens & Nelson|(|2007[) as 




( 6 ) 


where /?,„ and Rout are the disk inner and outer radii and 2^ and 
ec are the fluid element surface density and eccentricity respec¬ 
tively. Note Cc is calculated assuming the cell is orbiting the 
binary barycentre alone (a two body problem). Indeed, models 
with h - 0.05 show a higher value of dd at times both early in 
evolution and quasi-steady state (QSS), particularly for the low 
viscosity run. Our final disk eccentricity values lie between 0.03 
and 0.08, the spread matching closely that found in PN13. How¬ 
ever, we observe much larger values of eccentricity oscillation 
amplitude at QSS, which varies between models from 7.5 x 10“^ 
to 1.3 X 10“^ as opposed to a consistent value of 5.0 x 10“^ found 
by PN13. This may be due to our smaller value of initial sur¬ 
face density. The increase in eccentricity with aspect ratio is also 
seen in the radial profiles of the disk’s eccentricity (Fig.[^). The 
models with a larger aspect ratio (models B and D) have a larger 
average eccentricity than those with a smaller aspect ratio (mod¬ 
els A and C) from the inner boundary until about 2.5 au. 

We can also calculate the disk mean longitude of periastron, 
cbd, from the mass weighted average of the cell longitudes, 0 ;^: 




(7) 


The value oscillates about the binary periastron which itself does 
not change due to the analytic prescription for the stellar orbits. 
In Figure]^ cjd is plot as a function of time for Kepler-16. On 
cross examination with the eccentricity time evolution (Fig.[^), 
it is clear that there is some relation between the eccentricity and 
periastron longitude since the maxima in the eccentricity oscilla¬ 
tions matches the periastron longitude maxima. We discuss this 
link further in Sectionj^ The position angle of the disk centre-of- 
mass is calculated in the inertial frame and shows the circulation 
of the disk increasing in frequency until the disk precession pe¬ 
riod, Pd, settles at around 2,500 Pab at QSS. 

The eccentric morphology of the disk is apparent through 
both 1.20 (Figure [^1 and Cd- Figure shows the radial velocity 
which reveals a number of spiral features. In both the h - 0.03 
and 0.05 models of Kepler-16, the velocity map shows m - 2 
spiral waves with a larger radial extent for the higher aspect ratio. 
To better assess the level of asymmetry and presence of modes in 
the disk, a Fourier analysis of the surface density is performed. 

The disk surface density can be decomposed into modes 
described by azimuthal mode numbers m and frequency mode 
numbers I since disk disturbances caused by the binary can have 
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Fig. 4. Time evolution of disk mean longitude of periastron and the posi¬ 
tion angle of the disk centre-of-mass for Kepler-16 under both non self- 
gravitating conditions (black solid) and self-gravitating (red dashed) 
conditions. 


both an angular and time dependency. The surface density distri¬ 
bution can then be written as ( |Nixon & Lubow|2015| l 


nr, 0,t)^ Yd Tj exp[i(m0 - (8) 

l=-oo m=0 

where n,m{r) is a complex function and Q* is the mean motion 
of the binary, 27r/P.4g. [Nixon & Lubow| ( [2015| ) And that eccen¬ 
tric binaries produce modes with I + m. In particular, in their 
retrograde circumbinary disk simulations, highly eccentric bina¬ 
ries {eh > 0.6) produce powerful I - -I, m - 2 modes. We 
remove the time dependence of the disk disturbances by since 
our analysis is concerned with the eccentric Kepler-16 system 
{ch = 0.16) and because we focus only on a comparative study of 
the azimuthal mode strengths between simulations. We do this 
by choosing to discuss modes only with / = 0. This is done by 
performing the Fourier analysis of the surface density averaged 
over the period of a single binary orbit at QSS. This mode anal¬ 
ysis allows for the identification of the mode propagation and 
dissipation properties with changing disk parameters. 

In Figure [^ the strengths of the m - \ and m - 2 modes 
relative to the axisymmetric component is plot as a function of 
orbital radius. The results confirm that for the high aspect ratio 
ih - 0.05) models B and D, the contribution from the m - 2 
spirals is more significant in the outer disk than that seen for h - 
0.03. 

3.1.2. Kepler-34 

Disks surrounding Kepler-34 are significantly more affected by 
the binary, with the QSS Cd ranging from 0.1 to 0.15, typically 
three times more eccentric than Kepler-16 (Figure]^). Since the 
mass-weighted eccentricities decrease slightly with time by the 
simulation end, it is not clear from Cd if the disks have reached 
a steady state by 16,000 Pab- However, looking at the instan¬ 
taneous surface density profiles for a Kepler-34 run during the 


13 



























































S. Lines et al.: Modelling circumbinary protoplanetary disks 



Fig. 5. Fourier analysis of the surface density of all non self-gravitating 
Kepler-16 circumbinary disks. The transform is done on the time- 
averaged surface density, over a single binary orbit at 15,000 Pab, to 
consider only 1 = 0. The strengths of both the m = \ and m = 2 modes 
are normalised against the axisymmetric (m = 0) contribution and plot 
as a function of radius in the disk. 



Fig. 6. Instantaneous surface density profiles for Model H (Kepler-34) 
at three times near quasi-steady-state separated hy 640 Pab- 


last few hundred binary orbits as shown in Figure ^ it appears 
as though a steady-state has been reached. Our Kep er-34 results 
agree roughly with PN13, who find that typically, for the systems 
mutually covered, a quasi steady-state is achieved by 10^ Pab- 
As with Kepler-16 the average eccentricity in the Kepler-34 
disks once they have reached QSS is larger for higher values 
of aspect ratio. The large difference in Cd (Ae^ = 0.05) between 
models F and G is shown in Figure]^ to originate from the disk 
beyond the truncation radius where the surface density is much 
higher, and is therefore expected since Cd is a mass-weighted 
value. From 2.0 au, model F has a small positive eccentricity 
offset until 3.0 au, but there is very little difference between all 
models interior to this location. There is a high level of agree¬ 
ment between aspect ratio and a-viscosity models in the surface 
density profiles too; our simulations see f, sharing a similar value 
of 2.0 au across all models. These results agree with PN13. 

3.2. Self-gravitating disks 

In models I - L we enable disk self-gravity (DSG) for both 
Kepler-16 and Kepler-34. The o'-viscosity is fixed at 10“^ but the 
aspect ratio is varied between 0.03 and 0.05 to adjust the Toomre 
value and allow for direct comparison with models A, B, E and 
F. For Kepler-16, as seen in Figure]^ during the first 2,000 bi¬ 
nary orbits the disk responds almost identically to that of its non 
self-gravitating counterpart. The evolution up to QSS is subtly 
different however. The frequency of the eccentricity oscillations 
increases with gravity enabled suggesting that the disk circula¬ 
tion frequency increases. Due to the absence of strongly defined 
eccentricity oscillations in the Kepler-34 simulations, it is pos¬ 
sible only to differentiate between self- and non self-gravitating 
models by the eccentricity magnitude, and not the oscillations. 
We find a slightly reduced disk eccentricity with DSG enabled 
for the larger aspect ratio simulation. In Figure |7] the time av¬ 


eraged surface density profiles also reveal that a self-gravitating 
disk in Kepler-16 has no effect on the final disk structure, and 
almost no effect for disks around Kepler-34. 

To test at what disk mass self-gravity become important, So 
is increased by 5, 10 and 25 times. Self-gravity importance is 
therefore tested for Q ^ I but Q > 1 at all times. We define 
standard density models (So = 1 x 10^“* code units) as being 
0.5 'Lmmsn and those scaled up by a factor of 5, 10 and 25 as 
2.5 'Smmsn (Qmm - 28), 5 'Lmmsn (Qmin - 14) and 12.5 'Lmmsn 
(Qmin = 5) respectively. 

In Figure the evolution of the mean disk eccentricity for 
each of the four varying densities is shown, dd remains, for the 
non self-gravitating runs, at a constant 0.08 for each surface den¬ 
sity. When self-gravity is enabled, with the exception of 12.5 
'^MMSN^ falls off slowly with time and the rate of this fall 
off increases with increasing surface density. This leads to disks 
with higher surface densities having a slightly (5-25%) smaller 
eccentricity at the end of the simulation. The evolution of the ec¬ 
centricity in the self-gravitating 12.5 I-mmsn simulation is more 
dynamic due to an eccentricity fall off before increasing again. 

At standard densities we find that the oscillation frequency 
increases and amplitude decreases with self-gravity enabled. In¬ 
creasing the density does not change the oscillation frequency 
without self-gravity enabled, but the period of the oscillations 
decreases from 2000 Pab (0.5 'Lmmsn) to 300 Pab (12.5 I-mmsn) 
when it is included. 

Another measure of the effect of self-gravity is to identify 
the strength of the modes present in the disk. A fourier analysis 
of the surface density, identical to that done for the non self- 
gravitating Kepler-16 runs, is done. The Fourier coefficients are 
determined and normalised against the axisymmetric (m = 0) 
component. The strength of these modes normalised to the ax¬ 
isymmetric part, Om/fFo, is shown for each surface density in 
Figure]^ Again, a comparison is made between simulations with 
and without self-gravity included. 
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R (au) 

Fig. 7. Comparison of Kepler-16 disk eccentricity evolution and time 
averaged surface density for 0.5 'Lmmsn density disks with (/, J, K & L) 
and without self-gravity (A, B, E & F). Self-gravitating runs are shown 
with dashed lines. 


There is a strong contribution from both the m - I and m — 2 
modes, with the m - I mode always 20% or larger in the in¬ 
ner disk. It is less clear as to how the inclusion of self-gravity 
changes the strength of these modes. When self-gravity is en¬ 
abled in the 0.5 1,mmsn, 2.5 1.mmsn and 5 1,mmsn surface den¬ 
sity runs, there are only subtle differences between the strength 
of the modes. At 12.5 '^mmsn^ due to oscillations in the strength 
of both modes in the self-gravitating disk, the profiles for both 
the self-gravitating m - I and m - 2 modes do not trace their 
non self-gravitating equivalents. 

3.3. Surface density profile gradient 

The gradient of the initial density profile is a value which of¬ 
ten varies between studies, justified by the lack of observational 
data. In simulation U the density gradient is made shallower by 
changing the surface density exponent a-^ from 1.5 to 0.5. In Fig¬ 
ure the radially varying surface density is shown alongside 
the time evolution of the mean disk eccentricity. Changing the 
exponent has little effect on the truncation of the inner disk edge 
since I.peak adopts the same value. However, there is a significant 
difference in the mean disk eccentricity at all times, with aj; = 
0.5 obtaining a much lower eccentricity throughout. 

3.4. Boundary conditions 

Previous work on circumbinary disks sees the implementa¬ 
tion of two boundary conditions for the inner edge; namely 
rigid/reflecting and open. Open boundary conditions are chosen 



Fig. 8. Evolution of the mean disk eccentricity for standard density self- 
gravitating Kepler-16 run J and increased density runs 0-T (dashed 
lines). The corresponding simulation with self-gravity disabled is shown 
by solid lines. 


to be the more realistic scenario since matter can flow through 
the boundary and accrete onto the star(s) as would occur in a real 
system. The difficulty in these cases however is setting a well de¬ 
fined value for Vr and v^. Rigid boundary conditions are chosen 
to remove the enhanced mass loss experienced through the inner 
disk from an eccentric disk. The rigid approximation is, however, 
unrealistic and is often modified to include routines that remove 
reflected waves at the boundaries. We explore the influence of 
the choice of boundary condition on the disk by comparing the 
surface density and eccentricity between identical simulations, 
but with different inner boundary conditions. 

We find that there are small differences in both the disk ec¬ 
centricity and structure between the two. In Kepler-16 the sur¬ 
face density profiles show that including an open boundary keeps 
the peak density at the same location but increases the mass inte¬ 
rior to 1.peak- This effect, which is shown in Figure [TT]is stronger 
in Kepler-34 where the density plateaus at arounoL 1 au to 1.5 
au before further increasing to a peak at 2 au. We present our ex¬ 
planation for this in the Section The disk eccentricity reveals 
less of a difference between the two boundary conditions. In the 
main part of the disk the eccentricity is the same regardless of 
boundary condition although the mass weighted eccentricity is 
higher in the cavity. 
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Fig. 9. Fourier analysis of the surface density of Kepler-16 simulations testing the relevance of self-gravity with increasing Eq. Black lines 
correspond to the Fourier coefficients of the m = 1 modes, normalised to the axisymmetric m = 0 component. Red lines show the m = 2 modes. 
Simulations where self-gravity is enabled are shown by dotted lines. Fourier analysis is done on the time-averaged surface density over the period 
of one binary orbit from 15,000 Pab- 


4. Discussion 

In this paper, we have shown the results of a series of hydrody- 
namical simulations that explore the response of a circumbinary 
gas disk to both the Kepler-16 and Kepler-34 stellar binary sys¬ 
tems. We vary a-viscosity and aspect ratio, as well as probing 
the effects of disk self-gravity, initial surface density prohle, and 
the inner boundary condition. 

The majority of our simulations have reached steady state 
(a consistent value of e^) by 16,000 Pab- The steady-state na¬ 
ture of the disk is also reflected in the stationary form of the 
instantaneous, azimuthally averaged surface density during the 
last few thousand binary orbits (Figure [^, showing the balance 
between the truncation effects from tidal torques on the inner 
edge of the disk and the viscous replenishment of mass to the in¬ 
ner disk. The formation of the central cavity, which is shown by 
the movement of 'Lpeak in Figure!^ causes the inner edge of the 
disk to move outwards. By takingffie cavity edge to be the loca¬ 
tion at which the positive density gradient begins, we show good 
agreement with the Artymowicz & Lubow ( 1994[ ) gap estimate 
Pgap ~ 2 . 5 ^ 2 /,. 


The results for the surface densities and eccentricities of our 
non self-gravitating disks agree with those found by PN13. Sub¬ 
tle differences such as the eccentricity magnitude (Figurej^, par¬ 
ticularly in Kepler-34, may be explained by a number of factors 
that include our reduced outer boundary (from 5 au to 4 au) and 
the difference in initial surface density profiles. We certainly re¬ 
trieve higher disk eccentricities than found by |Kley & Haghigtr] 
ipour|(20T4ll who attribute their lower values in comparison with 


PN13 to the choice of boundary condition. |Kley & Haghighipour| 
( |2014| l claim that the rigid boundary condition used by PN13 
leads to a higher disk eccentricity. We And that the choice of a 
rigid boundary condition leads to only the slightest increase in 
disk eccentricity over the open boundary scenario in Kepler-34, 
and leads to lower disk eccentricities for a < 1.5 au in Kepler-16. 
Therefore, we suggest that it might be more appropriate to link 
the discrepancy in disk eccentricity to the choice of stellar binary 
for which we have shown significantly affects the disk structure 
and dynamics. 

The shape and strength of the binaries potential invariably 
changes with its orbital properties: mass ratio, eccentricity and 
semi-major axis. Therefore, it is likely that a circumbinary gas 
disk is only weakly affected by the relatively low eccentricity of 
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Fig. 10. Comparison of surface density and eccentricity evolution in 
Kepler-16 runs between density profiles of form S(r) oc Eq'" with aj 
= 0.5 (blue dashed) and 1.5 (black solid) 



Fig. 11. Comparison of open (dashed) and closed (solid) boundary con¬ 
ditions in Kepler-16 (black) and Kepler-34 (blue) via surface density 
and eccentricity profiles. 


Kepler-38 (e = 0.10) ( Orosz et al.|2012b| ) compared to the sig¬ 
nificantly higher eccentricity of Kepler-34 (e = 0.53). We would 
therefore expect the gas disk around Kepler-34 to adopt a higher 
eccentricity than a similar disk around Kepler-38. This theory is 
supported by the signihcantly larger disk eccentricities in sim¬ 
ulations around the eccentric Kepler-34 compared to Kepler-16 
where the stellar binary has a much smaller eccentricity. 

Increasing disk eccentricity with increasing stellar binary ec¬ 
centricity is a trend also seen in A^-body results. The planetesimal 
response in the Wbody scenario is dependent upon the forcing 
of stellar binary on the disk, and perturbation theory provides 
a value of eccentricity of which planetesimals tend towards. 


This forced eccentricity, Cf, is given by (Moriwaki & Nakagawa 
[2004l l 


e/ 


4 


Mb ah 1 + 3e^/4 




R^‘'\+3ell2' 


(9) 


where and Mb are the primary and secondary stellar masses, 
M, is the total binary mass, R is the distance from the binary 
barycentre and a* and Cb are the binary semi-major axis and ec¬ 
centricity, respectively. While it is clear that e/ varies with e* 
from Equation]^ how Cd is precisely affected by e* is not known 
since Kepler-34 and Kepler-16 differ in mass ratio as well as 
stellar binary eccentricity. The remaining question is do the gas 
and planetesimal disks adopt the same eccentricity; are e/ and 
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Cd equal? On average we hnd at 1 au, Cd ~ 0.2 for Kepler-16 (as 
per Figure ^), whereas e/ ~ 0.02. Kepler-34 has an even larger 
difference between its value of Cd ~ 0.45 and e/ ~ 0.002. 

There is a noticeable trend in Cd with aspect ratio but no 
clear change in the disk eccentricity or structure with varying 
cr-viscosity. We hnd considerably higher eccentricities in disks 
with a higher aspect ratio. This is explained by the increased 
level of communication in the disk huid due to the larger value 
of the sound speed, Cs- The higher Cj allows waves in the disk 
to propagate over longer radial distances before being damped. 
This results in an increased pitch angle of the spiral arms and 
increases the radial extent to which the waves can travel. Waves 
which unfold further out can deposit their energy over a larger 
percentage of the disk which increases the mass-weighted mean 
disk eccentricity, Cd (see Fig. |^. This effect is seen in Figure]^ 
as the enhanced propagation from a larger aspect ratio leads to 
a more significant contribution of m = 2 modes in the outer re¬ 
gions of the disk. The deposition of energy as a result of these 
unfurling spirals in the large aspect ratio runs, is also indirectly 
observed through the enhanced eccentric m - I modes in the 
outer disk, as the binary transfers energy to the fluid elements. 
Additionally, the radial eccentricity profiles of Kepler-16 show 
that a high aspect ratio leads to eccentricity increases in the outer 
regions of the disk. This can be seen in Figure [^, where the ec¬ 
centricity beyond 1.5 au in Kepler-16 is close to 0 for the low 
aspect ratio runs A and C, but remains greater than 0 for both 
runs B and D. This is also the case for Kepler-34 which has e = 
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Fig. 12. Comparison of surface density and eccentricity profiles 
in Kepler-16 runs between self-gravitating (dashed) and non self- 
gravitating (solid) for 0.5 Smmsjv (black) and 12.5 I.mmsn (red). The 
12.5 'Lmmsn surface density profile is normalised to the initial half 
minimum-mass solar nehula surface density. 


0.05 at 3 au for both high aspect ratio runs F and H, but close to 
0 for the low aspect ratio runs E and G. 

There is a correlation between oscillations in the disk eccen¬ 
tricity and the rotational dynamics of the disk from cross exam¬ 
ining the longitude of periastron in Fig ure and the disk ec¬ 
centricity in Figure |Pierens & Nelson ( 2013| l suggest that the 
amplitude of the oscillations is linked to pressure effects. This 
is reflected in our results by the increasing amplitude of eccen¬ 
tricity oscillations with larger aspect ratio. For a disk around an 
eccentric binary |Lubow & Artymowicz| ( |2000| l find the direct 
forcing from the m - I binary potential (Oij) on the disk can 
become important and the growth of the disk eccentricity is de¬ 
scribed (their Equation 2) as 


ed = -yUi)(l - 2fih)(abladf(l - e]) ’^D.asin'nj, (10) 

lo 

where Cb, fdb and a* are the binary eccentricity, mass parameter 
(MglMf) and semi-major axis respectively, is the Keplerian 
velocity of the disk and cr is the longitude of periastron of the 
disk relative to the binary. We do not retrieve the same magni¬ 
tude of eccentricity growth as predicted by Equation [1^ but this 
is likely due to the difference between the model’s dependency 
on ballistic particles in contrast to our fluid. The eccentricity 
growth and damping is well phased though with m, mapping the 
eccentricity oscillations in Kepler-16. Since the binary is placed 
on a fixed orbit cJj s m. Therefore by reading the evolution of 
the longitude of periastron of run B in Eigure during periods 
where zrr < 0, Equation gives ed > 0 and hence eccentricity 
growth. These periods are matched by an eccentricity increase 
for run B in Eigure]^. Similarly regions where vj > 0 and hence 
ed <0, a damping of eccentricity is seen in the eccentricity plot. 
The comparison of the phasing of ed with m for Kepler-34 is dif¬ 


ficult since the eccentricity oscillations in disks around this bi¬ 
nary are less pronounced. This validates the model further how¬ 
ever, as = 0 for Kepler-34, since the stars have equal mass 
(///, = 0.5). A similar effect is observed in the A-body case where 
the mass ratio of Kepler-34 leads to the removal of secular forc¬ 
ing (ef = 0 and dynamical forcing eff + 0) (|Paardekooper et al. 

[2OT2I 1. 

There is a discrepancy between the centre-of-mass position 
angle showing perfect circulation of the disk and the longitude of 
periastron showing a combination of imperfect circulation and li- 
bration. This may be explained by the mass weighted value of the 
longitude of periastron becoming contaminated by outer regions 
of the disk where higher order disturbances can cause the value 
of u to flip such that the orientation of the eccentric disk can vary 
between the inner and outer regions. If U) remained on the same 
hemisphere of the disk for all radii then wj would extend fully 
from -Ti to H-TT to show full circulation as per the centre-of-mass. 

We find a very minimal contribution from self-gravity to the 
dynamics of the disks. At standard densities, in both Kepler-16 
and Kepler-34 the eccentricity is marginally lower with gravity 
enabled, an effect observed by Marzari et al. ( |2009| l in their study 
of circumstellar disks perturbed by an s-type binary (circumpri- 
mary) configuration. Additionally, the oscillation frequency in ed 
increases with gravity enabled due to the increase in precession 
rate of the disk as seen in Eigure]^ Both these effects have only 
a minor impact on the disk’s density profile and its eccentricity. 
Therefore including disk self-gravity in simulations of disks at 
half minimum-mass solar nebula size is not necessary in mod¬ 
elling the disk and retrieving the correct eccentricity and struc¬ 
ture. 

We also find that self-gravity has no real influence on the 
evolution of the magnitude of the disk eccentricity or steady- 
state surface density at 2.5 and 5 times the minimum-mass solar 
nebula. At 12.5 times, self-gravity starts to become a necessary 
addition to the simulations to resolve the correct morphology of 
the disk. Enabling it at this high disk mass changes the struc¬ 
ture of the disk, as can be seen in the surface density plot of 
Eigure [1^ A series of small density humps in the outer disk is 
observed that occur at regions of increased disk eccentricity. A 
noticeable difference in the disk evolution when self-gravity is 
included is the change in frequency of the eccentricity oscilla¬ 
tions. The trend in final mass weighted mean disk eccentricity 
is to decrease with increasing surface density, with the excep¬ 
tion of the highest density mn. The disk eccentricity is probed 
further in Eigure where the eccentricity of the inner disk is 
greatly reduced with self-gravity enabled for 12.5 ^mmsn- The 
eccentricity then increases over the non self-gravitating run for R 
> 2.0 au. This may explain why the mass weighted value is larger 
at the end of the simulation, since it takes time to raise the ec¬ 
centricity of the outer disk. Including self-gravity therefore acts 
to slightly damp the overall eccentricity of the disk. Global self- 
gravitating modes can exist in the disk if precession due to the 
companion star can be ignored ( |Papaloiz^|2002[ ). However, it 
is likely that disk precession due to a massive binary companion 
cannot be ignored, in which case the global modes play no role 
of importance. 

The Eourier analysis performed to identify differences in the 
strength of the low frequency modes did not reveal any signifi¬ 
cant trends with increasing surface density, with or without self¬ 
gravity included (Eigure]^. With the exception of 12.5 'Lmmsn 
where including self-gravity is required to accurately model the 
m = 1 modes in the outer disk, it does not appear that including 
self-gravity is important when considering disks with Q > 14. 
Regardless of the surface density magnitude, the results reveal 
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the importance of including up to the m - 2 contribution of the 
surface density in fully describing the structure of the disk since 
outside the cavity, within the disk, the strength of the m - 2 
modes is often 20% of the background axisymmetric surface 
density. 

Our results indicate that the choice of boundary condition 
does not reflect strongly in the dynamical evolution of disks 
around both binary types. In both binary configurations, allowing 
material to flow through an open boundary increases the mass 
interior to the density peak. This effect is stronger in Kepler- 
34 than Kepler-16. Explanations for this include advanced wave 
damping routines operating in the rigid boundary case that could 
change the structure of the inner disk, or that the open bound¬ 
ary allows for a more realistic truncation timescale. There is 
little difference in the disk eccentricities beyond the truncation 
radius. These results are in contrast with |Kley & Haghighipour| 
( 2014| l who find much larger disk eccentricities when using a 
rigid boundary. This is probably to be expected since their shal¬ 
lower initial surface density profile (E oc and omission of 

a gap function to model the expected inner cavity means that a 
large amount of disk mass is located just exterior to the inner 
boundary. This means a large disk mass is perturbed close to the 
binary, increasing disk eccentricity. In this case, an open bound¬ 
ary is quick to remove any highly eccentric material close to the 
boundary itself. 

The final result from this work is the effect the initial sur¬ 
face density gradient, a^, has on the eccentricity. By changing 
aj; from 1.5 to a shallower gradient of 0.5, we retrieve a much 
smaller eccentricity at all times. Using - 0.5 reduces the fi¬ 
nal disk eccentricity by a half. This result is not surprising since 
the calculation of the mean disk eccentricity is mass weighted 
and therefore biased towards areas of high mass. For - 0.5, 
more mass is in the outer disk, as seen in the surface density pro¬ 
files of Figure Since the outer disk is less perturbed by the 
binary and has a lower eccentricity, disks with a shallower den¬ 
sity profile will have a lower mass weighted eccentricity. This 
explains why, alongside the varying binary parameters, Kley &| 
Haghighipour (|2014|l And considerably smaller values for their 


disk eccentricities than found by us and PN13. 


5. Summary and further work 

In this paper we have performed a number of hydrodynamical 
simulations of circumbinary gas disks that explore the binary 
configuration, disk aspect ratio, a-viscosity, inner boundary con¬ 
dition, initial surface density profile and inclusion of self-gravity 
on the evolution and quasi steady-state of such disks. 

We found that the choice of stellar binary has a significant 
impact on the eccentricity of the disk with Kepler-34 pump¬ 
ing up the mean disk eccentricity to three times the value seen 
in some simulations of Kepler-16. The truncation radius is also 
strongly affected by the binary configuration, with the more ec¬ 
centric Kepler-34 forcing the inner edge out to an average radius 
of 2.0 au, compared to 1.0 au in Kepler-16. There is no signif¬ 
icant correlation between the a-viscosity and the final surface 
density but there is a noticeable trend for the disk eccentricity to 
increase with increasing aspect-ratio. The eccentricity of the disk 
is strongly linked to the initial surface density gradient, which 
must therefore be taken into account when comparing work of a 
similar nature. At least for our initial surface density profile gra¬ 
dient value of aj; = 1.5, the choice of inner boundary condition 
did not have a substantial effect on the evolution of the disk for 
both Kepler-34, and Kepler-16, but may play a role in resolving 
the disk interior to the density peak around the truncation radius. 
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We did not see a strong influence of self-gravity in disks near 
minimum-mass solar nebula size, and saw only a minimal de¬ 
crease in mean disk eccentricity when self-gravity is included. 
Fourier analysis of the surface density to reveal the presence of 
modes also suggested that enabling self-gravity, even in disks 
where the surface density causes the Toomre parameter Q to ap¬ 
proach 1, does not change the structure of the disks. We saw a 
subtle difference however in the m — I mode strength and vari¬ 
ation in our high mass (12.5 ’^mmsn) model when self-gravity is 
included, suggesting that work on disk masses larger than this 
value should include self-gravitation in order to accurately re¬ 
solve the structure of the disk. This contrasts results showing the 
relevance of self-gravity in the circumprimary case since we did 
not And that it is a necessary inclusion in circumbinary disks that 
lie above the Toomre stability limit. 

The next goal is to take surface density and velocity profiles 
of a steady-state disk for both Kepler-16 and Kepler-34 for use 
in gas potential feedback on our A-body simulations of planetes- 
imal evolution and growth. From this work we found that many 
of our parameter choices have little influence on the Anal prod¬ 
uct, with the exception of aspect-ratio and initial surface density 
gradient which both have a significant impact on disk eccentric¬ 
ity. This will be an important aspect to consider in our future 
work since the eccentricity of the gas disk is a measure of the 
level of asymmetry in the gas disk potential. 

Future work on the hydrodynamical evolution of circumbi¬ 
nary disks will include the thermal response of the gas by includ¬ 
ing an energy equation. Additionally, since the evolution of the 
disks is sensitive to binary configuration, further investigations 
should look at changing the binary parameter space; eccentric¬ 
ity, semi-major axis and mass ratio. This may help to identify the 
link between dd and ef. 
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